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The Gran Sasso muon puzzle 
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We carry out a time-series analysis of the combined data from three experiments measuring the 
cosmic muon flux at the Gran Sasso laboratory, at a depth of 3800 m.w.e. These data, taken by the 
MACRO, LVD and Borexino experiments, span a period of over 20 years, and correspond to muons 
with a threshold energy, at sea level, of around 1.3 TeV. We compare the best-fit period and phase 
of the full muon data set with the combined DAMA/Nal and DAMA/LIBRA data, which spans the 
same time period, as a test of the hypothesis that the cosmic ray muon flux is responsible for the 
annual modulation detected by DAMA. We flnd in the muon data a large-amplitude fluctuation with 
a period of around one year, and a phase that is incompatible with that of the DAMA modulation at 
5.2o". Aside from this annual variation, the muon data also contains a further significant modulation 
with a period between 10 and 11 years and a power well above the 99.9% C.L threshold for noise, 
whose phase corresponds well with the solar cycle: a surprising observation for such high energy 
muons. We do not see this same period in the stratospheric temperature data. 



I. INTRODUCTION 

There has been a recent resurgence of interest in al- 
ternative explanations for the annual modulation sig- 
nal detected by DAMA/LIBRA, the dark matter di- 
rect detection experiment located at the Gran Sasso Na- 
tional Laboratory (LNGS), Italy. The DAMA/LIBRA 
data lU E], together with those of its previous incar- 
nation, DAMA/Nal [311], show a clear modulation (at 
8.9cr) that is consistent with the dark matter hypothe- 
sis both in period and phase. However, there has been 
widespread skepticism for the interpretation of this signal 
as evidence for dark matter direct detection. Alternative 
interpretations have been proposed: one of these is back- 
ground induced by cosmic muons, the flux of which also 
modulates annually with a peak in the summer, in the 
northern hemisphere, due to temperature fluctuations in 
the stratosphere. A (convoluted) mechanism by which 
the modulating cosmic muons might give rise to a sig- 
nal in the DAMA detector, involving intermediate spal- 
lation neutrons, was proposed in [Jj . This hypothesis has 
not been independently tested, and recently several ar- 
guments against it were put forward by DAMA [B] . Here 
we will examine more closely one of these arguments, 
namely the compatibility between the period and phase 
of the DAMA and muon annual modulations. 

Independent assessments of the compatibility of the 
DAMA signal with the cosmic muon flux, the latter 
taken from the LVD experiment [7j at LNGS, whose pe- 
riod of data taking coincided with the first 5 runs of 
DAMA/LIBRA, were carried out in [8] and [9], with con- 
tradictory conclusions. However, LVD is not the only ex- 
periment measuring the cosmic muon flux at the LNGS 
site; including the data from MACRO TU] and Borex- 
ino gives a 20-year span of muon modulation data, 
fully encompassing the time-span of both DAMA/Nal 
and DAMA/LIBRA. We analyse for the first time the 
combined data set and find an annual modulation whose 
phase is rather incompatible with DAMA's. Intriguingly 
for such high energy muons, we also see significant power 



at a period of just over 10 years, with a phase that rep- 
resents a close anticorrelation with the solar cycle. We 
present our results in Sec. [H] and discuss their implica- 
tions in Sec. IIIII 



II. RESULTS 

A. Muons and DAMA 

We begin by subtracting, from the data of LVD and 
Borexino, the average muon flux reported over the course 
of the experiment, in order to normalize them to the 
same baseline (MACRO already presented its data in this 
form). We then carry out a simple chi-squared fit of the 
combined data to a cosine of unknown amplitude, pe- 
riod and phase, marginalising over an added constant for 
each individual experiment, to allow for the effect of sys- 
tematic flux mis- measurements, as well as their different 
sensitivities for through-going muons^. The best-flt co- 
sine has a period of 365.9±0.2 (solar) days and a phase 
of 177.4±2.2 days (with respect to January 1st 1991). 
While these numbers are generally in good agreement 
with the fits carried out by the individual experiments, 
the goodness of fit, quantified by a value of chi-squared 
per degree of freedom of 7587/4244, is rather poor.^ This 
is unsurprising: while we expect the annual modulation 
of muons, which is directly related to temperature fluctu- 
ations in the stratosphere, to be periodic, this periodicity 



^ Neither the period nor phase of the leading behaviour change 
significantly on inclusion of these constants. 

^ Our assessment of the best fit parameters and uncertainties for 
each individual experiment are in good agreement with the values 
quoted by the collaborations themselves, with the exception of 
LVD's uncertainties in the period and phase which are an 
order of magnitude larger than ours, at fifteen days. Such a 
large shift, particularly in the period, would certainly be visible 
by eye over an 8- year time-span, and we find no evidence for this 
in the data. 
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is unlikely to be sinusoidal. Similarly we carry out a chi- 
squared fit to DAMA/Nal and DAMA/LIBRA data and 
find a period and phase in good agreement with those 
quoted by the collaboration themselves. We plot confi- 
dence limit contours for the muon and DAMA data sets 
in the 2D period-phase plane in Figj2] below. As pointed 
out in 9 , because the periodicity is allowed to vary, the 
size and shape of the contours are affected by the choice 
of time origin. We have verified, however, that the rela- 
tive (dis-)agreement between the two sets of data in chi- 
squared units is independent of the choice of the origin, 
as it should be, and is also relatively stable over the en- 
tire timespan of the experiments. We use the parameter 
goodness of fit [12], to quantify the level of compatibil- 
ity between DAMA and the muon data when measuring 
the period and phase. The p-value for the two data sets 
measuring the same parameters is 2.2 x 10~^, which cor- 
responds to a 5.2(7 tension between them. This supports 
the conclusion in Chang et. al. IQj of no strong correla- 
tion between the annual oscillation of cosmic muons and 
the DAMA signal, with two caveats. The first is that the 
mechanism by which the muons generate a signal in the 
DAMA detector does not significantly smear the phase of 
the modulation^ . The second is that we test for a correla- 
tion with the assumption of sinusoidal behaviour for the 
cosmic muon modulation, which is a rather poor one, as 
can be seen from the value of the chi-squared per degree 
of freedom. Indeed, when fixing the period to one year 
and extracting the phase for each independent year of the 
two datasets to test for the stability over time of the re- 
sults, we found that the yearly DAMA results are always 
compatible within their uncertainties with the average 
value over the length of the experiment. On the other 
hand, the cosmic muon data shows a larger yearly dis- 
persion evidencing again the poorness of the sinusoidal 
approximation for the muon data. However, we expect 
that a more sophisticated statistical analysis of the data 
sets, that does not rely on this assumption, will yield a 
similar conclusion."* 



B. Muons and the solar cycle 

Even more interesting perhaps are the results from 
our systematic tests for the existence of a subleading, 
longer-term variation in the muon flux, such as that no- 
ticed by Blum |5] in the LVD data. We fit the combined 
muon data to a sum of two cosines, with unrelated and 



^ A simple smearing by a Poisson distribution to account for the 
stochastic nature of this process, proposed in 8 , is insufficient 
to overcome the tension between the two data sets, for smearing 
by any distribution of reasonable width. 

* Note that even taking LVDs quoted uncertainties at face value 
does not significantly change the degree of discrepancy between 
phases of the cosmic muon flux and DAMA data since MACRO 
and Borexino data are enough to provide a A.Ta tension. 



unknown amplitudes, periods and phases, marginalising 
over a constant shift parameter for each experiment as 
before (see Fig. [l] for best-fit curve) . We find the lead- 
ing annual modulation almost unchanged, but the chi- 
squared for the fit improves by over 250 units on addition 
of the second cosine, with an amplitude of 0.40±0.03%; 
a period of 10.7±0.3 years and a phase of 1880±50 days 
(corresponding to a first maximum in March 1996). As a 
quantitative measure of the significance of the subleading 
periodicity in muon data, we plot a Lomb-Scargle peri- 
odogram [T5HT5] . for the full cosmic muon data set. Like 
an inverse Fourier transform, this prescription separates 
a periodic signal into its harmonic constituents, but is 
tailored to work even for unevenly-spaced data. It has 
well-understood statistical properties, for example it is 
known to yield an exponentially-falling power when the 
input is random Gaussian noise. Given the large dispar- 
ity between the size of the error bars in the data from 
the different experiments, however, we find it necessary 
to weight the data points by their individual error bars, 
as proposed in [16j . thus making the information in the 
periodogram more analogous to that obtainable using a 
chi-squared fit [Uj. We use an oversampled set of fre- 
quencies based on the natural frequencies of the data 
set in order to obtain a good resolution on the result- 
ing periodogram, and estimate the threshold power at 
which noise can be excluded at the 99.9% C.L. using the 
prescription provided in [17 . We generate 10,000 sam- 
ples of random Gaussian noise, with identical spacing to 
the muon data, and with the same variance and error 
bars. We then compute the weighted Lomb-Scargle pe- 
riodogram for each sample, for our standard set of over- 
sampled natural frequencies, and use the distribution of 
the power of the highest peak in each plot as a measure 
of the "false-alarm" probability. Our results can be seen 
in Fig. |3ja) . The dominant feature is a peak at 1 year 
(truncated due to the range of the plot). There seem 
to be many subdominant peaks, including one around 
10 years, although it is unclear how to interpret them, 
since it is only possible to come to a statistically rigorous 
conclusion about the dominant peak in a given L-S peri- 
odogram. We therefore repeat the above procedure after 
filtering out the annual modulation, by subtracting from 
the data a cosine with the best-fit annual parameters. 
The resulting periodogram can be seen in Fig. [Sjjb) , with 
the limit for exclusion of noise at the 99.9% confidence 
level (now given the presence of an annual modulation 
with the pre-defined properties) shown as a dotted line. 
In the subtracted periodogram one can clearly see in the 
cosmic muon flux a harmonic component with a period of 
10.4±0.3 years^. Many of the small-period spikes present 
in the periodogram for the full data set have vanished af- 



^ This is in fact the same period that is obtained using a chi- 
squared fit, in the absence of the constant term, added to account 
for systematics. 
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FIG. 1: Combined cosmic muon data from MACRO, LVD and Borexino, with monthly binning for clarity, after subtraction of 
the mean measured flux at each experiment, as well as of an additional constant, determined for each individual experiment by 
the chi-squared fit. This constant, which accounts for systematic differences in the experiments' sensitivities to cosmic muons, 
has a negligible effect on the best fit period, phase and amplitude for the annual modulation. The best fit to the sum of two 
independent cosines yields the solid line in the figure (periods, phases and amplitudes detailed in the text). 
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FIG. 2: Confidence-limit contours for period and phase of 
best-fit cosine for cosmic muon flux from MACRO, LVD and 
Borexino, and DAMA data. The straight black-dashed lines 
delineate the expected period and phase for a dark matter 
signal. See text for note concerning changes in time origin. 



ter harmonic filtering of the annual modulation, leading 
us to conclude that they were aliasing peaks, or some 
other artifact of the irregular data spacing or sampling 
frequencies used. The subleading periodicity in the muon 
data, as seen in both the chi-squared fit and the L-S pe- 
riodogram, displays a strong correlation with the solar 
cycle. We also plot for comparison the weighted Lomb- 
Scargle periodogram for the monthly-averaged sunspot 
data in the same period, taken from |18| . and originally 
derived from data by the Solar Influences Data Analy- 
sis Center in Belgium in Fig. Islc). We see a dominant 



peak corresponding to a period of 12.6±0.1 years, and 
what looks like higher harmonics of this fundamental fre- 
quency. Caution must be exercised in interpreting the 
fitted periods and uncertainties presented in this subsec- 
tion: the solar cycle is known to have a rather variable 
period, making the cosine fit an inadequate description 
of the data, as reflected in the large chi-squared values. 
This does not, however, preclude the use of these fit val- 
ues to compare two data sets under the hypothesis of a 
correlation between them. 

In order to test for a possible correlation between the 
secondary modulation in the Gran Sasso cosmic muon 
data and the sunspot data, we again make use of the pa- 
rameter goodness of flt, with parameters extracted from 
a chi-squared fit of both the muon and sunspot data to 
cosine functions with a relative phase of tt. While the 
fitted phases are in agreement, we find there is a 4.7(7 
tension between the fitted periods, which is not very en- 
couraging. Notice, however, that in the flts, the sunspot 
data are weighted by their variance, which in the limit of 
low statistics is strongly dependent on their absolute val- 
ues, making data taken during minimums of solar activ- 
ity dominate the fit. Because of this, the fitted period is 
driven to large values by the unusually long and deep so- 
lar minimum around 2008. By contrast, the correspond- 
ing muon data was taken mostly by Borexino, which has 
comparatively large error bars. Thus, muon data give 
more weight to earlier parts of the solar cycle, which fit 
better with smaller periods. Indeed, redoing the fit af- 
ter rescaling the sunspot error bars such that the relative 
size of the error bars (and hence their weights) are the 
same as in the muon data reduces the tension between the 
two data sets to 2.\a. Furthermore, as mentioned above, 
a sinusoid with a constant period is a particularly bad 
model for sunspot activity (522/242 for the chi-squared 
per dof), which is cyclical rather than periodic. Perform- 
ing a fit to the sunspot data using a cosine with a period 
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(a) Weighted Lomb-Scargle periodogram 

for cosmic muon flux variation. The 
dominant feature is the peak at 1 year. 
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(b) Weighted Lomb-Scargle periodogram 
for muon flux residuals after filtering out 
annual variation. 
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(c) Weighted Lomb-Scargle periodogram 
for sunspot data. 



FIG. 3: Weighted Lomb-Scargle data showing subleading 10.4-year period in cosmic muon flux variation, and correlation with 
sunspot data, for the period January 1991 to April 2011. 



that varies linearly in time instead results in a very ac- 
curate description of solar data over the two cycles in 
question, with an improvement in the chi-squared per 
dot (to 249/241), and a fitted "period" which varies from 
8 to 13 years. With this phenomenologically-motivated 
fitting function, we find that the discrepancy between 
both data sets reduces to less than la, seemingly imply- 
ing a close correlation between them.^ This correlation 
is rather puzzling given that the cosmic muon flux is be- 
ing measured at 3800 m.w.e. below the earth's surface, 
which corresponds to a threshold muon energy of 1.3 TeV 
at sea level. 



The variation of the cosmic ray flux with the solar cycle 
is a known phenomenon, one that is now understood to 
be due to larger magnetic fields, and increased turbulence 
from the solar wind at a solar maximum, deflecting low 
energy cosmics, thus preventing them from reaching the 
earth. According to the common lore [iniUnjj however, 
one should not expect to see this effect persist at ener- 
gies larger than tens of GeV for the primary cosmic, and 
hence also for its decay products. The persistence of this 
effect to muons of energies larger than O(TeV) is likely 
not due to long-term changes in the stratospheric tem- 
perature, since we do not see the corresponding period 
in the effective temperature periodogram (computed as 
detailed in [21], using atmospheric temperature records 
taken at the nearby Pratica di Mare station, from the 
Integrated Global Radiosonde Archive [12 )■ There is 
possibly a more subtle mechanism at work; coming up 
with a suitable candidate, however, will require a detailed 
analysis of the interplay between cosmic ray propagation 
and atmospheric effects, and is beyond the scope of this 
paper. 



® In contrast, we were unable to flnd an alternative function that 
provided a better fit to the annual modulation of muons. 



III. DISCUSSION 

We combined the measurements of the cosmic muon 
flux from three Gran Sasso-based experiments, MACRO, 
LVD, and Borexino, and analysed the resulting 20 years 
of data in the light of claims that cosmic muons might 
somehow be responsible for the 8.9cr annual modulation 
measured at DAMA. In fltting the muon and DAMA data 
to a sinusoidal variation of unknown amplitude, period 
and phase, we flnd the two data sets have periods that 
are compatible, but their phases are in conflict at 5.2a. 

Without a working model by which muons might fake a 
signal in the DAMA detector, a possible role for cosmic 
muons in DAMA's annual modulation cannot be ruled 
out by these arguments alone. However, it seems chal- 
lenging to flnd a mechanism that can bridge a phase dis- 
crepancy of 5.2a between the two datasets. Alternatively, 
DAMA might be measuring the resultant effect of two 
annual modulations with slightly different phases, one of 
which results from contamination by cosmic muons. 

We subsequently separated the muon data into its har- 
monic components using a Lomb-Scargle periodogram. 
We found in addition to the annual modulation, a sub- 
dominant modulation of period just over 10 years with a 
power well above the 99.9% C.L. for noise, and a phase 
that is anticorrelated with the solar cycle. This result 
was confirmed using a chi-squared fit. 

A correlation between such high energy muons and the 
solar cycle goes against the common lore: one might ex- 
pect energetic cosmics that produce these muons to be 
unaffected by the presence of the larger solar magnetic 
fields and stronger solar winds of a solar maximum. This 
puzzling observation is unlikely to be due to contami- 
nation by cosmic neutrinos, whose flux at the depth of 
the Gran Sasso lab (3800 m.w.e.) would be subdomi- 
nant to the muon flux, and too small to account for this 
effect [53]. Moreover, we found no evidence for such long- 
term modulation in the effective stratospheric tempera- 
ture close to Gran Sasso. The reason for the persistence 
of this effect to high-energy muons possibly stems from 
some complex interplay of atmospheric effects and sec- 
ondary cosmic production, but the detailed modelling of 
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these effects is beyond the scope of this work. 

Independent tests of a potential correlation between 
high energy cosmic rays and the solar cycle should al- 
ready be possible at a number of long-running exper- 
iments, including underground detectors such as Super- 
Kamiokande, and terrestrial experiments such as AGASA 
or HiRes. Additional facilities exist that are currently 
taking cosmic ray data, like IceCube and MINOS for 
muons, or Extensive Air Shower experiments such as the 



Tibet Air Shower Array and Argo-YBJ. In a number of 
years, these will have collected enough data to probe the 
relevant time scales, allowing us to explore the depen- 
dence of any modulation with the depth, latitude and 
longitude at which the observations were recorded, as 
well as its energy and flavour-dependence. Finally, satel- 
lite experiments can directly probe the primary cosmic 
rays, yielding crucial information for our ultimate under- 
standing of this effect. 
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